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ABSTRACT 

Recent work in galaxy formation has enlightened the important role of baryon physics, to solve the main problems 
encountered by the standard theory at the galactic scale, such as the galaxy stellar mass functions, or the missing 
satellites problem. The present work aims at investigating in particular the role of the cold and dense molecular phase, 
which could play a role of gas reservoir in the outer galaxy discs, with low star formation efficiency. Through TreeSPH 
simulations, implementing the cooling to low temperatures, and the inclusion of the molecular hydrogen component, 
several feedback efficiencies are studied, and results on the gas morphology and star formation are obtained. It is shown 
that molecular hydrogen allows some slow star formation to occur in the outer parts of the discs. This dense and 
quiescent phase might be a way to store a significant fraction of dark baryons, in a relatively long time-scale, in the 
complete baryonic cycle, connecting the galaxy discs to hot gaseous haloes and to the cosmic filaments. 
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1. Introduction 



Galaxies: evolution — Galaxies: ISM — Galaxies: spiral — Galaxies: star formation 



Numerical simulations of galaxy formation have now 
reached a high degree of sophistication, and include increas- 
ingly detailed physics, from star formation rates treated 
with sub-grid recipes based on the Kennicutt-Schmidt law 
(KS), to more direct processes triggered by Jean s instabil- 



ities ([Stinson et al. 2006 Hopkins et al. 2011), feedback 



that is treated either through supernovae heating, or ki- tures (iBlanchard et al 
netic impulse, momentum- driven flows d ue t o stellar winds 
dSales et al!]|2010[ [Ostriker fc^ hetty 2QTTt, and in some 
cases, multiphase gas ( |Maio e t al. 2007; tinedin et al.|2009| ). 
However, big unknowns remain as free parameters in the 
process: first the nature of dark matter, and its behaviour, 
leading to well- identified problems at galaxy scales, such 
as predicted cuspy profiles (while only cores are observed), 
or predicting a large number of dark satellites around each 
Milky- Way type galaxies. Second, a large unknown is also 
the nature and location of the missing baryons: observation- 
ally only less than half of the baryons have been identified 
(e.g. [Fukugita fc Peebles 2004 J , and this missing baryonic 
component is certainly gaseous, given th e results of mi- 
crolensing (e.g. Wyrzykowski et al. 2009[ ). This gas might 
be hot (of the order of millions of degrees) or cold, and 
most of it must be located in the intergalactic space, not 
to overpredict rotation curves. 

A significant fraction of dark baryons should also exist 
in galaxies, and could reside under the form of cold gas, 
dense enou gh to be in the molecular phase (e.g. Pfenniger^ 
Combes|| 1994; Gr enier et al.||2005| |Bournaud et al.||2007| ' 
Langer et al.||2010J . A reservoir of cold gas may help to 
moderate star formation, and explain almost stationary 



star formation his tories in lat e- type galaxy discs as ob- 
served today (e.g. Wyse 2009). When the thermal evolu- 
tion of gas is studied during galaxy formation, an immedi- 
ate conclusion is that almost all the baryons should have 
cooled and condensed in dark matter potential wells, where 
they are supposed to form stars, which are not seen. This 
overcooling problem can be solved through re-heating the 
gas, e ven before its inflow in the dark matter bo und struc- 

However 



Dave et al. 



2001) 



the star formation efficiency might be lower in globally 
low gas density environments, such as the outer parts of 
galaxies, and cold gas reservoirs are another solution to ex- 
plore. The time-scale spent by the gas in the cold phase is 
largely unknown, and this could change the baryon cycle. 
Cosmological simulations are unable to derive the abun- 
dance of cold gas, by lack of spatial resolution, leading to 
reduce density dynamics, underestimating the cooling lo- 
cally. 

In the recent years, the cold mode gas accretion to as- 
semble mass onto galaxies has been revived versus the hot 
mode gas accretion. It has been realised that only part of 
the gas was heated to the virial temperature of a structure 



and inflowing 


^, and most of the gas 


(Keres et al.| 


2005|. The cold gas 


for small ha] 


oes, lower than 10^^ 



^0, 



is inflowing along the cosmic filaments, while for massive 
systems, the hot quasi-spherical mode dominates. There is 
also a redshift dependence on the importance of the two 
modes, the cold mode being more dominant a t hig h red- 
shift {z large r than 2). |Birnboim fc Dekel (2003) and Dekel 
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fc Birnboim| ([2006) found similar conclusions, using one- 
dimensional simulations and analytic arguments. The limit 
between dominating cold and hot modes corresponds to a 
baryonic mass close to the scale separating the red sequence 
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and the blue cloud in gala xy populations (Kauffmann et al. 
2003 Baldry et al. 2004), and the origin of bimodality in 
observed galaxies has been connected to these two differ- 
ent theoretical assembly modes. Even in massive structures, 
cosmological simulations have shown that cold filaments of 
gas may penetrate their hot gaseous halo (Birnboim et al. 
2QQ7| [Dekererar||2Q09| ). 

A further well-known problem in the baryon cycle and 
galaxy formation is to reproduce correctly the galaxy stel- 
lar mass function. If the large amount of predicted dwarf 
galaxies are suppressed by supernovae feedback, assumed 
to expel gas out of their small potential wells, gas is still 
accreted in massive systems, and forms constantly grow- 
ing, blue star-forming galaxies at the massive end of the 
mass function. Solutions have been proposed through AGN 
feedback (either quasar mode through heating and out- 
flows, or radio mode through powerful radio jets), with 
limited success however. In particular, the gas expulsion 
is more efficient again in low- mass galaxies. Besides, even 
when ejected with a speed larger than the escape veloc- 
ity, gas is re- accreted onto galaxies, being braked not only 
through gravity but also by hydrodynamical interactions 
with the halo gas. The re-accretion occurs on time-scales 
of the order of one Gyr, and ha s been called halo fountain 
by Oppenheimer & Dave (2008). According to the environ- 
ment, gas outflowing from low or intermediate mass haloes 
can be re-accreted by more massive companion galaxies, 
a pr ocess which has b een dubbed the intergalactic foun- 
tain ( Keres et al.||2009 ). Globally, the frequent re- accretion 
of the ejected gas limits the efficiency of the AGN feed- 
back, maintaining some overcooling problem, and keeping 
massive galaxies active in star formation until late times, 
with no definite quenching. It is possible to introduce in 
the baryon cycle a cold gas reservoir, with a sufficiently 
long time-scale in this phase to account for the observed 
properties of galaxies. In the present paper, we want to 
explore this possibility, beginning by isolated galaxy simu- 
lations, and considering in a companion paper galaxies in 
their environment. We are taking into account a cold dense 
molecular phase in the gas component, which can cool at 
a temperature much lower than the lO^K usually assumed 
for "cold gas" in cosmological simulations. According to the 
baryon physics adopted, and in particular the amount of 
feedback, we will examine the stability of the galaxy discs, 
and the efficiency of star formation, in order to derive the 
order of magnitude of the time-scale spent by the gas in 
this cold component reservoir. 

Thanks to the always increasing computing capacities, 
there has been in recent years a growing number of simula- 
tions taking into account the cold phase. These works come 
from two different domains, with different goals: 



structure. These simulations are in general multiphase, 
with the goal to answer questions about the cloud for- 
mation, the star formation efficiency The galaxy models 
are often not self-consistent, with no transfer of angu- 



1. 



2. 



Cosmologically oriented work, tending to higher spa- 
tial resolution, of the order of 100 pc; their main in- 
terest is to determine whether the small scale structure 
impacts the large scale tranfer of angular momentum, 
the stability of discs, formation of bulges, concentration 
of the galaxy, transformation of galaxies from late to 
early type, whether feed back ca n destroy dark matter 
cusps etc. dMaio et al.|20 07; G nedin et al.|2009||Schaye 
eFar]|2QTQl |Murante et al.pOlOllBournaud et al.||2010| 
Hopkins et al.||201lD; 

Star formation oriented work, tending to enlarge the 
field of view from molecular clouds to the whole galaxy 



dSlyz et al. 20051 Tasker & Bryan 


2006, 2008 Dobbs 


& Bonneh 120081 IDobbs et al.||2011 


Shetty & Ostriker 


20081 Wada et al.||2011). 



Some works are at the transition between the two scales, 
and when the simulations are unable to resolve the cloud 
fragmentation scale, have tried to represent the multi- 
phases of the gas using sticky particles, moving ballistically 
in the potential, and able to collide and coagulate to form 
a whole mass spectrum of clouds (e.g. ,Semelin fc Combes] 
2002[ [Booth et al.||2007| [Revaz et al.]|2Q09p . 

In the present work, our goal is to deal with the multi- 
phase of the gas in a self-consistent way, allowing the gas 
to cool down to 100 K, the average temperature of the qui- 
escent dense molecular phase, and to follow its interaction 
with the other phases through cooling/heating, star for- 
mation and feedback, through a Tree-SPH code. The H2 
molecule formation will be simplified by a recipe related to 
the density, metallicity and self-shielding of the gas. The 
star formation recipe is based on the Kennicutt-Schmidt 
law, and core-collapse supernovae feedback is taken into 
account. 

The numerical techniques and the initial conditions of 
the simulations carried out are described in ^ together 
with the details of the baryonic physics adopted to deal 
with cold molecular hydrogen, star formation and feedback, 
^describes the results of the simulations. The influence of 
the feedback and inclusion of molecular hydrogen cooling 
is detailed. Our conclusions are drawn in Q 



2. Numerical techniques 



We use the Gadget-2 code (Springel 2005) that computes 
gravitational and hydrodynamical forces. Gravity is com- 
puted by a Tree algorithm with a Barnes Hut opening 
angle of 0.7, taking only the monopole moments of the 
grav itational f ield in to account for computational reasons 
(see |Springel (|2005[) for details). We take the same con- 
stant softening length e for all particle types (stars, gas 
and dark matter) . Hydrodynamics is treated with a Smooth 
Particle Hydrodynamics (SPH) algorithm, with an individ- 
ual smoothing length h computed such that a constant mass 
is contained in a sphere of radius /i, and that is allowed to 
be as low as 0.1 e. 

We add baryonic physics to this code, as described in 
this section. 



2.1. Cooling 

Gas in the interstellar medium (ISM) loses internal energy 
by radiation due to physical processes occurring at micro- 
scopic scales. This cooling can be encapsulated in a cooling 
function representing the corresponding volume rate of en- 
ergy loss, and depending on the chemical composition of the 
gas, its density and temperature. In simulations of galaxies, 
cooling is often considered only down to a temperature of 
~ 10^ K, as the most efficient cooling processes are due to 
ionised hydrogen and helium, elements that are almost only 
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2.1.2. Metal-line cooling below 10"^ K 




Fig. 1. Cooling functions nor malised by the atomic hydro- 
gen number density, following S utherland fc Dopi ta (1993) 
for high temperatures (T > 10"^ K, and Maio et al. (2007) 
for low temperatures (T < 10^ K), as a function of temper- 
ature and metallicity Z = [Fe/H] (Z = is solar). In the 
code, we use more exact density-dependent cooling func- 
tions for the low temperature part, cf Sec. 2.1. 



Following Maio et al. ( 2007 ), we compute cooling functions 
of a few metals below lO^K: CII, 01, Sill and Fell, as they 
are the most abundant heavy elements released by stars in 
the ISM and are thus the main ingredients for cooling, con- 
sidering their collisions with H atoms and electrons. For 
"low" densities obtained in galaxy simulations, the pop- 
ulations of energy levels needed to compute the cooling 
rates differ from the Boltzmann populations of the Local 
Thermal Equilibrium (LTE). The populations here depend 
on the abundances of species that can collide with the metal 
and allow it to change of energy level, which makes the rate 
of energy loss per metal depend on these abundances (con- 
trarily to the LTE case for which this rate depends only 
on temperature). The resolution of the equations governing 
the populatio ns is included in th e code, using the quantum 
data given inlMaio et al.| (|2007|), and the cooling functions 
are computed using the obtained populations and the tran- 
sition probabilities, as in Mai o et al.^ (2007) . 

We choose the same abundances as in [Sutherland fcl 
Dopita ( 1993[ ) (solar abundances and primordial ratios) for 
the metals we consider here and a low electronic fraction of 
— lO"^)^ assuming the cold gas is quasi neutral. 

The cooling curves for uh = 1 cm~^ and different metal- 
licities Z = [Fe/H] are plotted in Figure [l] 



present in the atomic form below 10^ K. This usually im- 
plies the coldest dense gas is at an equilibrium temperature 
of about 10^ K in these simulations, as dense gas experi- 
ences heating pressure forces and shocks if it reaches lower 
temperatures, and the gas is prevented from collapsing fur- 
ther. We include cooling by metals, molecular hydrogen and 
HD down to 100 K, a more realistic temperature floor for 
the ISM. 

In dense media, cooling by metals might dominate at 
solar metallicity. However, metals are abundant only in the 
central parts of galaxies, and radial abundance gradients 
are observed in giant spiral galaxies, with an exponential 
decline of about 0.6 dex in 10 kpc o n average (e.g. Henry 



Howard 



1995| [van Zee et aL||1998D . Therefore the outer 

parts of galaxies are reminiscent of the primordial galaxies, 
with low metal abundance. Cooling through the H2 and HD 
molecules is then important, and could change the physics 
of star formation in these regions. Ultraviolet observations 
by GALEX have shown that star formation can be active 
at large radii, much farther than the optical radius R25, 
in regions where Hq, observations are not able to reveal 
moderate-age populations of stars. UV-bright discs extend 
up to 3 — 4 times the optical radius in about 30 % of spira l 
galaxies ( |Thilker et ar||2005[ [011 de Paz et ar]|2005[ |2QQ7 ). 
The presence of molecular hydrogen and star formation in 
outermost discs of spirals is of prime importance to study 
cold gas accretion, w hich is considered one key factor i n 
galaxy evolution (e.g. |Keres et al.||2005[ |Dekel et aL||2009D . 



2.1.1. H, He and metals cooling above 10^ K 



We use the c ooling functions computed by [Sutherland fc 
Dopita| ( 1993) for a plasma in collisional ionisation equilib- 
rium above 10^ K. These cooling functions include cooling 
due to ionised H, He and metals and are tabulated for dif- 
ferent metallicities. 



2.1.3. H2 cooling 

We are interested in studying the influence of molecular hy- 
drogen cooling on galaxies, especially in regions where met- 
als are not abundant. We tabulate the LTE cooling function 
as a function of temperature from the Boltzmann energy 
levels populations and transition probabilities. For low den- 
sities for which the cooling function is lower than the LTE 
cooling function, we use data from Glover fc Abel| (2008) 
for a 3:1 ortho-para ratio of H2. We assume H2 cooling be- 
low 10^ K takes place in an almost neutral medium and 
thus consider collisions with the following species: other H2 
molecules, H atoms and He atoms. Collisions with H atoms 
are often the only one taken into account, however, as we 
are interested in describing media with a high fraction of 
H2, taking into account the collisions with H2 and He is nec- 
essary to compute a cooling that would otherwise vanish in 
regions poor in atomic hydrogen. 

The final expression for the volume cooling rate (in 
ergcm"^ s~^) we use in the code is: 



^H2 



A 



H2,LTE 



1 



A 



li2, LTE 



(1) 



nHAHo,H+nHo Ae 



? +^HeAH2 ,He 



Ah2,lte is the LTE cooling rate per H2 molecule, in 
ergs~^, and Ah2,h, Ah2,H2 AH2,He are the cooling func- 
tions in ergcm^s"^ due to collisions with respectively H, 
H2 and He, for low densities . This formula, similar to what 
Hollenbach fc McKee (1979) used, interpolates between the 
low density and high density regimes: it makes the influ- 
ence of collisional processes increase with respect to the 
LTE cooling as the density decreases. We plot in Figure |2] 
this cooling rate for a hydrogen nuclei number density of 
1 cm~^, and different mass ratios of molecular hydrogen 
over the total hydrogen component. 
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We do not consider cooling due to collisions of H2 
molecules with metals as collisional coefficients are not well 
determined. 




Fig. 2. Cooling rates of H2 for a hydrogen nuclei number 
density of 1 cm ~^ as a function of temperature for different 
H2 to total hydrogen mass ratios xu^ . For high fractions of 
H2 and at densities for which the cooling rate depends on 
collisions, the cooling can become slightly more important 
just below 10^ K as the mass fraction is reduced, because 
the contribution of collisions with atomic hydrogen domi- 
nates in this domain of temperatures. 



2.1.4. HD cooling 

We also consider HD cooling. Despite its lower abundance 
than molecular hydrogen, the molecule HD can help cool 
the gas because of its permanent electric dipole. The dipole 
transitions probabilities are orders of magnitude greater 
than the quadrupole ones similar to the H2 quadrupole 
transitions, which makes the cooling per molecule signif- 
icantly greater. We used quantum data to compute the 
LTE cooling function in a similar way tha n for H2. Fo r 



low densities, we took the function of Lipovka et al.| ( 2QQ5[ ). 
This cooling function takes into account the collisions of HD 
with H atoms and only dipolar transitions, which is justi- 
fied by their much bigger contribution to the cooling than 
quadrupole ones. We followed Glover & Abel (2008) by as- 
suming the cooling function is simply proportional to the H 
atoms density for low densities, so that the cooling function 
per molecule is, for low densities, nHAHD,H(^H = lcm~^) 
We write the volume cooling rate, as for H2, as: 



A 



HD 



^HD- 



A 



HD,LTE 



Ah 



(2) 



T'HAHD,H(nH = l cm-3) 



Ahd,lte is the LTE cooling rate per HD molecule, in 
erg s~^, and AHD,His the cooling function in erg cm^ s~^ due 
to collisions with H atoms. 

In the simulations, we assume that nno = lO^^nRo- 



2.1.5. Cooling algorithm 

The interstellar gas component is modelled as an ideal gas 

5 

with an adiabatic index 7 = -. Gadget-2 integrates the 

o 

P 

evolution of the specific entropy A = ^ rather than the in- 
ternal energ y, for computational cost and accuracy reasons 
developed in Springel & Hernquist (2002). The evolution of 
the specific entropy of a particle i is governed by: 



dt 



1 



1 ^ — 



(3) 



where A is the volume cooling rate in erg cm^ s ^ (the 
sum of the contributions described in above), Hij is a fac- 
tor accounting for the artificial viscosity and Wij is the 
symmetrised SPH kernel. 

As the cooling time can be lower than the time-steps 

equal to the minimum of yjj^ (where r] is an accuracy 

parameter we keep at 0.025 as in the default Gadget-2 
parametrisation, e is the gravitational smoothing length, 
and a is the acceleration), and the Courant time limita- 
tion for gas particles, we use an implicit cooling scheme 
to stabilise the resolution. We solve iteratively the implicit 
equation: 



Ai{t^ dt) = Ai{t)^ 



7-1 



A(p,(t),T,(t+ dt)) 



P^(t)^ 
1 ^ 

-^m,H,,v,,.V,iy,, I dt (4) 



where dt is the considered time-step. 

The relation between temperature and specific entropy 



is: 



,7-1 



(5) 



T^iXi — - with Xi the number fraction of the 
rrip 

species i is the mean particle weight. For a gas composed 
of atomic and molecular hydrogen and helium, neglecting 
the metals contribution: 



where ji 



3-4 



(6) 



X 



^Hnuci is number density of hydrogen nuclei and X 
is the hydrogen nuclei mass fraction that we set to 0.76. We 
do not compute the ionisation fraction — in the simu- 

^^nucl 

lations. We assume the gas is quasi neutral below 10^ K. 
Above that temperature, we assume there is no molecular 
hydrogen and the gas is fully ionised (H+ and He^+). The 
mean particle weight being higher for a neutral gas than for 
an ionised one, we set to 10^ K the temperature of the gas 
in the specific entropy range for which the temperature as- 
suming an ionised gas is below 10^ K but the temperature 
assuming a neutral gas is above it. 
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Table 1. Resolution 



e moM m^ 
[pc] [Mo] [Mo] [Mo] 
"Too sTfW IAW 2.5 10^ 



Table 2. Resolution: gravitational softening and particle 
masses. 



2.2. Resolution 

Distinguishing physical effects from numerical ones can 
be difficult if the resolution of a simulation is not well 
adapted. The number of particles setting the mass reso- 
lution, the gravitational softening length and the hydrody- 
namics smoothing length need to be chosen so that they 
do not introduce unwanted artefacts while allowing for the 
simulation to run, and at a reasonable pace. 

We want to describe the gas down to low temperatures, 
but if the temperature is too low, the Jeans length and 
Jeans mass of the gas will be more poorly resolved, which 
can lead to numerical a r tefact s ( ,Bate fc Burkert^^l997 ). 
Schaye fc Dalla Vecchia (2008) use an effective equation 
of state at high densities that makes the Jeans mass in- 

2011) take a 



dependent of density, while Hopkins et al. 
density-depend ent pressure floor to ensu re the Jeans length 



is resolved, and Bournaud et al. (2010) use an equation of 
state that mimics the effect of cooling at all temperatures 
and a temperature floor at high densities. We take a tem- 
perature floor of 100 K for all densities. We note that since 
we simulate disc galaxies, the gas is stabilised by the rota- 
tion of the disc. 

The number of particles should ideally be the largest 
possible to have a good mass resolution, which allows the 
Jeans mass to be well resolved. We run simulations with 
1 200 000 initial particles: a third are gas particles, a third 
are stellar particles and a third are dark matter particles. 
Table [2] shows the particle masses we have in the particular 
case of our galaxy model that will be detailed in |2.6[ New 
stellar particles have a mass mnew * = J^^ depending on 



the parameter Ng^^ described in 2.4, the number of stellar 



particles produced out of one gas particle. 

The gravitational softening length e depends on this 
number of particles: a too small value for a given number 
of particles introduces unphysical two-body relaxation in 
media that are collisionless (stellar components of galaxies 
and dark matter haloes) , while a too large one decreases the 
spatial resolution by "blurring" density features and does 
not allow the Jeans length to be gravitationally resolved. 
Gravitational softening lengths e are usually taken as scal- 
ing with the mean inter-particle distance, therefore as the 
number of particles to the power | for a 3-dimensional sim- 
ulation. We take a softening length e = 100 pc for all parti- 
cle types. This particular value is derived from the GalMer 



par 

simulations (eg'Di Matteo et al.|2007 ) which took ten times 
less particles for a softening of 280 pc. The GalMer simu- 
lations were isothermal, at 10^ K, and since here the tem- 
perature reaches down to 100 K, we take a softening length 
that is a little inferior to the cubic root, while still allow- 
ing for an efficient computation on a few tens of computing 
cores. 

The value of the gas smoothing length h should be 
small to allow for a good density resolution. Gadget-2 uses 



variable smoothing lengths: the densities and smoothing 
lengths are computed together, so that the mass contained 
in a sphere of radius is a constant. We take a minimum 
h equal to a tenth of the gravitational softening. This min- 
imum can have subtle consequences on the structure of the 
gas. In the presence of strong cooling and with no sufficient 
star formation or efficient feedback, a number of particles 
can be "stuck" at this minimum, making the correspond- 
ing regions increasingly denser with no possibility of getting 
more diffuse. This can be computationally very demanding 
and lead to large domain work-load imbalance. 

2.3. H2 fraction 

Our goal is not to determine the molecular abundance 
through a detailed chemical scheme, which would be too so- 
phisticated for our present approximated treatm ent (which 
does not t a ke int o account radiative transfer). [Krumholz 
Gnedin ( 20111 ) have compared different methods, and 
proposed a simple analytic approximation to compute the 
mass fraction fn^ (ratio of the molecular hydrogen mass on 
the total hydrogen nuclei mass). We follow this recipe in 
writing: 



/h2 



4 1-h0.25s 



with: 



s = 



ln(l + 0.6x + 0.01x^) 



(7) 



(8) 



Tc is the dust optical length and x is a scaled UV flux 
over the hydrogen nuclei number density. is set to zero 
when s > 2. 

The dust optical length is rc = • The dust cross 

section per H nucleus ad is set, using the reference Milky 
Way value, to (Jd = Z'10~^^cm~^, with Z' the metallic- 
ity we normalise by the solar metallicity. /in is the mean 
mass per H nucleus. H is a column density H = ob- 
tained from a local scale height: L = ~^^~y We compute 

{Vp)i^ the gradient of the density of the particle i by: 
Piiyp)i = Y.j'^j{pj - pi)ViW(rij,hi). The scale height 
takes the variation of the density into account: it increases 
with density but is inversely proportional to its gradient, so 
it is lower in the case of large gradients encountered on the 
outer parts of density features like clumps or spiral arms. 
It is thus well adapted to compute the column density used 
to determine the shielding from radiation. 

The X is computed using the new stars radiation field 
G. We assume stars whose age is inferior to 10 Myr ra- 
diate in the Lyman- Werner band and can dissociate the 
II2 molecules, and do not consider that older stars radiate 
in this band, so they don't have any effect on the frac- 

G 

tion of II2 in our simulations. We take x ^ , with 

^Hnuclei 

a scaling factor depending on the resolution of the simu- 
lations, and tuned so that the obtained II2 fractions are 
consistent with observations. The radiation flux produced 
by stars decreases with the distance r to a star in ^, as 
the gravitational field. We insert the computation of the 
flux received from stars younger than 10 Myr by one gas 
particle in the Gadget-2 gravitational tree functions, as it 
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is similar to computing the gravitational force due to the 
young stars only and can be easily done by adding a range 
of variables representing this contribution only. 

2.4. Star formation 

We implement star form ation in a s tochastic way that re- 
produces a Schmidt law (Katz][l992 ): 



dp* 
dt 



dt 



Pg 



(9) 



where is the volume stellar density, pg the volume 
gas density and tf[ is the free fall time: tff = ^2Gp 

is the star formation efficiency per free fall time. 

The Schmidt law can be enforced by giving a gas particle 
a probability to spawn a star at each time-step At: 



(l-e-^^*) 



(10) 



This implementation means that a gas particle of mass 
VTig can spawn a star particle of mass with this proba- 
bility at each time- step, the mass of the gas particle being 
then reduced by the amount of m*, until there is no more 
mass in the gas particle. The number of Ng^^ of star par- 
ticles created by gas particle is a compromise between a 
good mass and time resolution of star formation and CPU 
cost: if several stellar particles are created from a single 
gas particle, star formation will be smoother, temporally 
better resolved but the total number of particles can in- 
crease significantly, which slows down the code. In the case 
where Ng^^ is greater than one, we note that this spawning 
scheme implies that gas particles will have different masses, 

-.The 



either the initial one m^, or a smaller multiple of 
density and smoothing length are computed by Gadget-2 
so that the mass contained in a sphere of radius hi is fixed, 
hi being the smoothing length of the particle i, and we were 
careful about modifying the algorithm to take into account 
the different masses, so that this mass condition is still sat- 
isfied. 

We use common selection rules for the particles that are 
allowed to spawn stars: they must have a density higher 
than a threshold nHmin^ must have a temperature lower 
than a maximum temperature T^ax and must be in a con- 
verging flow (V.v < 0). We set nnmin = 10~ ^cm^. The 
threshold can be increased to allow for an ISM with more 
density structures, but if it is too high, the mass resolution 
must be increased for the Jeans mass to still be resolved. 
We use Tmax = 30 000 K, which happens not to be selective 
as only a few diffuse particles can reach these temperatures 
in our simulations. 

2.5. Feedback 

Our simulations include core-collapse supernovae kinetic 
feedback. We consider a supernova explosion releases the 
canonical value Esn — lO^^erg. A fraction of this en- 
ergy is given to the ISM, while the rest is radiated away. 
Each new stellar particle of mass inputs an energy 
-^input = Q^fb^sN^^ where esN is the mass fraction of formed 
supernovae in the new stellar particle multiplied by the 
canonical energy £^SN = lO^^erg. We consider this mass 



fraction is 0.01, which is of the order of the fractions ob- 
tained from commonly used initial mass functions (such as 
a Salpeter IMF with slope —1.35 and lower and upper limits 
of 0.1 and 40 Mq, with stars more massive than 8 M© con- 
sidered to be supernovae), so we have esN = lO^^erg/M©. 

Each neighbour i of a new star particle receives an 
energy weighted by its distance to the new stellar particle: 



Engbk^(kkoUo)' 



^input 



(11) 



SO that the sum of the energies given to the neighbours 
is £^input- If the newly created stellar particle has left a rem- 
nant gas particle, no feedback energy is given to the rem- 
nant. We recompute the smoothing length and neighbours 
list at the position of a new stellar particle, considering only 
the neighbouring gas that has not been turned into stars 
(and not considering a possible gas particle remnant at the 
exact position of the new stellar particle). Therefore the 
sum of the fractions involving kernels is indeed unity, and 
the mass of the gas affected by each supernova explosion is 
a constant. 

The neighbours are given a velocity kick directed 

along the line joining the stellar particle and the neighbour, 
and away from the new stellar particle. 

The number Ng^^ of stellar particles created out of one 
gas particle has an influence on the distribution of feedback 
energy. Feedback is more gradually input for a larger Ng^^. 

2.6. Initial Conditions 

We consider the case of a giant Sb galaxy. 

The gaseous and stellar discs follow Miyamoto-Nagai 
density profiles. The density of the gas disc is, in cylindrical 
coordinates: 



agR^ + (ag + 3^^^T7^) (ag + ^ ^ h^)^ 



5 

and the stellar disc density is: 



(12) 



Pd{R,z) = 



47r 
- (ad -+ 



■3^ 



hi) (a, + Vi^)' 



R^ + (ad + y^z^+f^)y {z^+hjr 



(13) 



The Miyamoto-Nagai profiles are chosen because the 
associated potential is analytic and the initial velocity dis- 
persions can thus be easily computed, but relaxation makes 
the profiles quickly exponential as observed. 

We use the gravitation routines of Gadget-2 to obtain 
the forces acting on particles. We compute the circular ve- 
locity of disc particles using these gravitational accelera- 
tions, and add an analytic asymmetric drift correction to 
have a more realistic velocity profile. 
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The radial velocity dispersion is derived from cr^(r) = 
3.36GQS(r) ^ ^^^.j^ q ^^iq Toomre parameter that we set to 

1 for both discs, S the surface density and k, the epicyclic 
frequency derived from the potential. The azimuthal ve- 



locity dispersion is obtained from 



with Q 



the angular speed, and the vertical dispersion is set by the 
isothermal equilibrium of the disc. 

The spherical stellar bulge and dark matter halo have 
Plummer density profiles. The bulge density is thus given 
by: 



Pb{r) 



3M5 
47rr5 



and the halo density is: 
3M^ 



Ph{r) = 



(14) 



(15) 



All the masses and characteristic lengths are specified 
in Table El 

The velocity dispersion of the spherical components is 
chosen to be isotropic and is derived from the second mo- 
ment of the Jeans equation. The radial dispersion is thus: 



CTl{r) 



1 

p{r) 



p(«)-dn 



(16) 



The velocity curve for this analytic model, with the con- 
tributions of the different components, is shown in Figure |3] 
The rotation is due mainly to the stellar components near 
the center of the galaxy, and to dark matter at large radii. 

We further prepare the initial conditions by letting the 
galaxy evolve for 300 Myr with only gravitational forces 
included for any type of particles. This allows us to start 
the simulations with dynamically relaxed discs that will 
especially exhibit no annuli instabilities. The initial surface 
densities of gas and stars are shown in Figure [4j and a 
snapshot of the gas is shown on Figure [5] 
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Fig. 3. Initial rotation curve. Top solid (blue) thick line: 
total rotation curve. Dashed (purple) line: bulge. Dotted 
(green) line: stellar disc. Dash-dot (green) hue: gas. Solid 
(red) line: dark matter halo 
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Fig. 4. Initial surface densities. Dotted line: gas. Dashed 
line: stars. 
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Fig. 5. Initial density map of the gas disc. The box size is 
[40 kpc X 40 kpc] 



Table 3. Sb galaxy parameters 



_Sb 

in Mr: 



1.7 10' 



4.5 10' 



Mb 
1.1 10' 



_JAg 

0.9 10' 



Th dd hd 



n 



da 



hn 



in kpc 12 5 0.5 1 11.8 0.2 



3. Simulations and results 

We assume there is a metallicity gradient in the gas. We 
take a central metallicity Z = [Fe/H] such that Z{R = 0) = 
Zq +0.5 and assume the metallicity decreases of 1 dex per 
10 kpc: 



Z{R) = Z{R = 0) - 



R[kpc] 
10 



(17) 



where R is the cylindrical radius. The influence of the 
cooling by metals will thus decrease with the distance to 
the centre of the galaxy. For simulations with no molecular 
hydrogen, we expect the outer regions to have fewer density 
features than the central ones and to be warmer. If however 
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we add some molecular hydrogen, depending on its fraction 
and on the details of star formation and feedback, the gas 
may be able to form clumps also in the outer regions and 
form more stars there. 

The gas has an initial temperature distribution peaking 
at a few hundreds of kelvins. The temperature floor we ap- 
ply is more exactly a specific energy floor, corresponding 
to a temperature of 100 K for purely atomic gas, and go- 
ing up to 186 K for gas with hydrogen present only in the 
molecular form. 

For all the simulations of the paper, the star forma- 
tion efficiency per free-fall time is set to — 0.1, and the 
number density threshold for star formation is nHmin = 
10~^ cm^. The simulations presented here all have a num- 
ber of stars spawned by gas particle Ng^^ = 4. 



3.1. Simulations without H2 

We first perform simulations with no molecular hydrogen, 
with varying feedback efficiencies afb- We have four dif- 
ferent feedback efficiencies: either no feedback (offb = 0), 
Qffb = 1%, afb = 10% or Qffb = 40%. 

Density maps of the gas discs for these runs are shown 
in Figure [6] at three simulation times: 0.5 Gyr, 1 Gyr and 
3 Gyr. The gas can reach smaller values of local volume den- 
sity when feedback is included, as star formation in dense 
regions inputs some energy in the ISM and prevents it from 
getting denser, slowing down the star formation. For the run 
without feedback, the central part of the galaxy exhibits a 
strong density contrast on the snapshots at 0.5 Gyr and 
1 Gyr, many small clumps and thin spiralling filaments can 
be seen. If feedback is switched on, we can see only a few 
clumps for a feedback efficiency afb = 1% and smoother spi- 
ral features, and no clumps for = 10% and afb = 40%. 
A central thin bar is formed in all the cases. After 3 Gyr, 
the bar is still clearly visible for the higher feedback effi- 
ciencies, but is less obvious otherwise, because the central 
parts have been depleted from gas by star formation. 

Figure [7| shows the specific energy- number density his- 
tograms after 0.5 Gyr of evolution, a time at which gaseous 
discs are actively forming stars in all the simulations. On 
the top and right of each plot, the marginal probability 
density functions (PDFs) of respectively hydrogen nuclei 
number density and specific energy are shown. The spe- 
cific energy PDFs show the fraction of gas at T :^ 10^ K 
(corresponding to the higher specific energy concentration) 
increases with feedback, while the cold dense gas fraction 
decreases. On the density PDFs we can see that gas reaches 
smaller maximum densities for higher feedback efficiencies 
and there is an increasingly high fraction of diffuse gas. In 
the low feedback runs, a significant fraction of the dense gas 
has a temperature close to the minimum: this is the dense 
gas of the central parts, subject to metal-line cooling. The 
fraction decreases for higher feedback efficiencies, because 
of the dissipation of energy by feedback, making the dens- 
est gas of the simulations warmer. If the gas is heated by 
pressure forces, viscous shocks or feedback, its temperature 
does not reach much beyond 10 because of the stronger 
H, He, and metals cooling it undergoes above. The diagonal 
branches observed on the left of each plot account for gas 
that, away from the centre of the disc, is subject to very 
little cooling because of the metallicity gradient, and thus 
cools down adiabatically. 




Fig. 7. Specific energy-density histogram of the gas after 
0.5 Gyr of evolution. The 2D histograms and marginal 
ID histograms are all mass weighted and normalised. The 
smaller specific energy value corresponds to a temperature 
of 100 K, and the top accumulation of gas is at 10^ K. From 
top to bottom: runs with an increasing feedback efficiency. 
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Fig. 6. Projections of the gas density after 0.5 Gyr, 1 Gyr and 3 Gyr of evolution. Box sizes are [30 kpc x 30 kpc] for 
face-on views and [30 kpc x 10 kpc] for edge-on views. Each row corresponds to a feedback efficiency indicated on th8 
left. T he column integrated density scale is the same for all plots. Done with the visualisation software SPLASH ( [Price 
20071). 
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We plot in Figure [S] the time evolution of the total mass 
of gas present in the simulations. All the gas is originally 
in the disc but can leave it under the effect of gravitational 
heating or stellar feedback. The characteristic time of con- 
sumption of the gas increases with feedback efficiency be- 
cause of the moderating effect of feedback on gas density, 
regulating star formation, and the curves have increasing 
horizontal asymptotes y-values. This is due to gas expelled 
from the disc, and also to the inability of the gas to reach 
the star formation threshold. The gas in the outer parts of 
the disc remains diffuse with almost no star formation in all 
cases, as the lower abundance of metals does not allow the 
gas to cool down enough to form stars. More details con- 
cerning the star formation will be presented in the following 
section when compared to simulations including molecular 
hydrogen. 



1.0 



0.8 



0.6 



!5 0.4 



0.2 



0.0, 





afb=40% 

— «fb=0% 









4 5 6 
Time [Gyr] 



Fig. 8. Time evolution of the mass of the gas component. 
Solid (blue) line: run without feedback. Dashed (green) 
line: run with a feedback efficiency afb=l%. Dash-dotted 
(red) line: run with afb=10%. Dotted (black) line: run with 
afb=40% 



since it depends on the physics at very small scale, below 
our resolution, and on complex radiation transfer, through 
gas clumps and associated dust. The global mass fraction 
first decreases with time because new stars are formed and 
contribute to the dissociating radiation field and becomes 
stable when the SFR becomes almost null (see the solid line 
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of Figure 15 for the evolution of the gas mass). Figure 
shows the surface density of H2 and atomic hydrogen gas 
HI after 0.5 Gyr of evolution for the case x x 50 that we 
choose. Such a surface density profile is simila r to some 
observed profi les of local galaxies d escribed in [Young fc; 
Scoville (1991), and more recently in Bigiel et al.r(2008). 




Fig. 9. Global H2 mass fraction versus time for afb 
The UV flux increases from top to bottom. 



3.2. Simulations with H2 

We now study the impact of the inclusion of H2 on the gas 
physical state, star formation and structure of the discs. 



3.2.1. H2 fraction 

We first set the feedback efficiency to afb = 1% and tune 
the X factor in the H2 fraction so that we find a global 
H2 fraction that is consistent with observations while hav- 
ing enough molecular hydrogen fraction to study its ef- 
fects. Our metallicity gradient implies that H2 will tend 
to be less present in the outer parts of the galaxy where 
the metallicity and dust abundance are low and the gas is 
less shielded from Lyman- Werner radiation. However, this 
is compensated by the lower star formation rate in these 
regions, which decreases the ambient Lyman- Werner lumi- 
nosity. We plot the evolution of the total mass fraction of 
H2 in Figure |9] for four different UV flux scaling factors. The 
exact exposition of the gas to UV flux is not well known. 




20 
K [kpc] 



Fig. 10. Radial distribution of the surface density of H2, 
atomic gas HI and total hydrogen gas after 0.5 Gyr for 
Qffb = 1% and the selected UV scaling factor. 
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Fig. 11. Projections of the gas density (top row) and H2 density (bottom row) after 0.5 Gyr, 1 Gyr and 3 Gyr of evolution. 
Box sizes are [30 kpc x 30 kpc] for face-on views and [30 kpc x 10 kpc] for edg e-on views. The column integrated density 
scale is the same for all plots. Done with the visualisation software SPLASH ( Price|[2QQ7 ). 




1.5 2.0 
Time [Gyr] 

Fig. 12. Fraction of cold g as as a function of time for two 
different feedback efficiencies Offb = 1% and Offb = 10%. 
Solid lines: with H2. Dashed lines: without H2. 



3.2.2. Gas physical state 

Figure 11 shows the aspect of the disc for = 1% at 
the same simulation epochs than in Figure [6) whose second 
row is the simulation for the same feedback efficiency but 
without H2. Cooling by collisions of atomic hydrogen with 
metals in the purely HI simulation or mainly by H2 in the 
now H2 dominated central part make this region similar 
in both simulations, the difference is in the outer parts of 
the galaxy for which metal-line cooling is poorly efficient 
because of our assumed metallicity gradient. In this case, 
the gas remains diffuse with no other cooling processes, 
but the inclusion of H2 allows the gas to be clumpier: we 
see density features that were absent in the purely atomic 
simulation. The surface density of H2 is also plotted. It 
indeed follows the density features of the gas: the clumps 
and filamentary structures can be seen in H2. 

The specific energy-number density histograms of 



Figure 13 and the marginal PDFs show higher fractions of 
gas in the cold dense phase than for the corresponding feed- 
back efficiencies without molecular hydrogen. There is no 
clear diagonal branch because there is some efficient cool- 
ing in the whole disc. The specific energy floor corresponds 
to gas between 100 K and 200 K (the exact minimum tem- 
perature depends on the molecular gas fraction), while the 
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Fig. 13. Specific energy-density histogram of the gas after 
0.5 Gyr of evolution for runs with H2. The 2D histograms 
and marginal ID histograms are all mass weighted and nor- 
malised. The smaller specific energy value corresponds to 
a temperature of 100 K for atomic gas, and the top accu- 
mulation of gas is atomic gas around 10^ K. From top to 
bottom: runs with an increasing feedback efficiency. 




Fig. 14. PDFs of hydrogen number density as a function of 
galactocentric radius after 0.5 Gyr of evolution. Top plot: 
run with no H2. Bottom plot: run with H2. Horizontal black 
line: threshold density for star formation. 



higher specific energy concentration of gas corresponds to 
gas at a temperature of about 10^ K. We separate the gas in 
two ranges of specific energy below and above 10 (km/s)^, 
which corresponds to a temperature between 1000 K and 
2000 K depending on the molecular fraction. We study the 
evolution of the fractions of cold and warm gas (gas lying 
below or above this threshold) depending on the feedback 
efficiency, and with or without the inclusion of cooling by 
molecular hydrogen. The majority of star formation hap- 
pens in the first Gyrs in all the simulations (especially when 
the feedback efficiency is low and stars form quickly) , so we 
focus on this period and plot the cold gas fraction as a func- 
tion of time on Figure [12] for two feedback efficiencies. All 
the gas is initially in the cold phase. For a given feedback 
efficiency, the cold gas phase represents a much higher frac- 
tion of the gas if H2 cooling is taken into account. Without 
H2, stars are formed from warmer and more diffuse gas, 
reducing the star formation efficiency. Feedback decreases 
the fraction of cold gas: the kinetic energy given to particles 
in dense star forming regions is transformed into thermal 
energy by pressure forces and viscosity, all the more as the 
feedback efficiency is high. 

As we are especially interested in studying the state of 
the ffas as a function of galactocentric radius, we plot the 
histograms of gas density versus radius in Figure 14 for the 
atomic and molecular simulations having afb = 1% , after 
0.5 Gyr of evolution. The effect on density is clear: the gas 
exhibits density peaks at larger radii and has a mean higher 
density with the inclusion of H2, with gas denser than the 
star formation threshold even at large radii. 



3.2.3. Star formation 

The effect of H2 cooling on the star formation efficiency is 
visible in Figure [15] and Figure [16] Figure [15] shows the total 
gas mass evolution with time for Offb = 1%. The depletion 
time of the gas decreases when molecular hydrogen cooling 
is included. 

The star formation efficiency is similar in the central 
regions, but the inclusion of H2 allows for star formation in 
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Fig. 15. Evolution of the gaseous mass of the galaxy. Solid 
(blue) line: run with H2. Dashed (green) line: run without 
H2 




15 20 25 

R [kpc] 

Fig. 16. Cumulative star formation rate as a function of 
galactocentric radius, averaged on the first Gyr. Solid lines: 
run with H2. Dashed lines: run without H2 



the outer parts too. Having inserted the formation time of 
stars in our simulation outputs, we can track star formation 
spatially. Our outputs are temporally spaced by 10 Myr. 
We define the SFR as being the mass of stars formed dur- 
ing 10 Myr divided by this time, which is similar to the star 
formation rates obtained from observations in H^. We plot 
the cumulative SFR averaged on the first Gyr as a func- 
tion of radius for these two simulations and other feedback 
efficiencies in Figure [T6j and we indeed see that about the 
same amount of star formation occurs in the central parts, 
but molecular hydrogen starts playing a role at large radii. 
The difference occurs at a larger radius for the simulation 
with no feedback, which can be explained by the already 
high clumping in central parts of the disc, making star for- 
mation very efficient even without H2. 

Figure [Tt] shows the projected density of stars formed 
since the beginning of the simulation for Offb = 1%: the disc 
of new stars is more extended if we include H2. There is 
a few stellar clumps, and also a very clear stellar bar that 



is maintained in the stellar component after a few Gyrs in 
both cases. Clumps of young stars can be seen on the den- 
sity maps: they follow the gas clumps. These stellar clumps 
gradually lose energy and are eventually absorbed by the 
central bar. 

We further study the star formation by drawing 
Kennicutt-Schmidt (KS) diagrams representing the surface 
density of SFR as a function of the gas surface density. As 
we are limited in mass resolution for star formation (stellar 
particles of a fixed mass of ~ 10^ Mq are created stochas- 
tically, completely differently from some smooth star for- 
mation), in order to have a significant amount of data to 
study, we add data points corresponding to 50 snapshots, 
from t=200 Myr to t=700 Myr. We use a polar grid with a 
given number n^^ of bins in cylindrical radius R and a given 
number n^i of bins in azimuthal angle 0. Using this kind of 
grid allows for a more uniform signal/noise repartition than 
with an orthogonal grid, and optimises the number of new 
stars per cell. 

On Figure [T8| we have plotted KS diagrams for simu- 
lations without feedback, and with different feedback effi- 



ciencies. These can be compared with |Agertz et al. (2012) 
plots of azimuthally averaged KS diagrams of disc galaxies 
for different feedback intensities. Very similarly, we observe 
a global diminution of the SFR surface density when we in- 
crease the feedback strength. The figure quantifies how, on 
average, for the same gas surface densities, the SFR is lower 
with higher feedback. The feedback makes the gas more 
diffuse and destroys clumps. Two cells containing the same 
amount of gas but different fractions of diffuse gas (cells 
are bigger than the clumps size), will have different star 
formation efficiency. This also explains the smaller scatter 
in the simulations with feedback: as the gas is more homo- 
geneous, the relation between surface densities of SFR and 
gas is better determined. In Figure [Tsj lines of constant gas 
depletion time are indicated. The gas depletion time is de- 



fined as tdep = 



^Gas 
^SFR 



. It can be seen that the high SFR and 



gas surface density regions of the galaxies have a depletion 
time as low as a few hundreds of Myrs in the simulations 
with no or little feedback, while the low SFR and density 
regions have depletion times up to 10 Gyr. The outer parts 
of the disc with a low surface density form stars much less 
efficiently than the central parts. 

We also see a difference for low surface densities between 
atomic and molecular simulations. The SFR is varying more 
linearly with molecular gas than with atomic gas. And the 
low surface density regions form stars more efficiently with 
molecular gas. This is explained by the fact that H2 cooling 
allows the gas to be locally denser. It especially allows it to 
be more concentrated in the disc plane as can be checked in 
the edge-on projections. Then the gas is denser in volumic 
density, and forms stars more efficiently, at a given surface 
density. 

For simulations including H2, the SFR is shown sepa- 
rately as a function of atomic and molecular components, 
on the bottom of Figure 18 The atomic hydrogen surface 
density is confined to low values for our galaxies and the 
SFR-HI diagrams show a large scatter because HI is too 
diffuse to track the star forming gas. The fact that H2 is a 
better tracer of star formation is also found in observations 
of nearby galaxies (e.g. Bigiel et al 2008). 
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Fig. 17. Projection of the density of stars formed during the simulations after 0.5 Gyr, 1 Gyr and 3 Gyr of evolution. 
Top row: feedback efficiency afb = 1% and no H2. Bottom row: feedback efficiency afb = 1% and H2. 



3.2.4. Vertical structure of the disc 




Time [Gyr] 

Fig. 20. Fraction of gas beyond 1 kpc from the disc plane. 
Solid lines are runs with H2 and dashed ones are runs with- 
out H2. Blue lines: no feedback. Green lines: afb = 10% , 
Red lines: Offb = 40% 



The inclusion of molecular hydrogen has a significant 
impact on the vertical structure of the disc, first because the 
cold and dense gas concentrates in the disc middle plane, 
and second because of the impact of gas clumping on the 
distribution of feedback energy. Christensen et al. (2012a,b) 




afb=40% with H2 

-- Qifj^ =40% without H2 
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Fig. 21. Vertical mass profile of the disc for radii R > 
15 kpc for runs with afb = 40%. Solid line: run with H2. 
Dashed line: run without H2. The two vertical lines mark 
the characteristic vertical heights zi/2- 



include non-equilibrium formation of H2 , self-shielding and 
dust shielding of both HI and H2 in galaxies extracted 
from cosmological simulations and explore the influence 
of including H2 formation, for a fixed feedback efficiency. 
Similarly to their results, we find that the introduction of 
H2 makes the outer parts of discs denser. But also it allows 
the gas to be more efficiently expelled from the plane, since 
the feedback has a stronger effect towards dense regions: in 
our feedback scheme, particles get velocity kicks weighted 
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Fig. 18. K-S diagrams. Left column: no feedback. Right column: Offb = 10%. Top row: simulations without H2. Second 
row: simulations with H2. Bottom plots: SFR versus the HI only and H2 only surface densities for simulations with H2. 
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Fig. 19. Vertical structure of the gas. First row: characteristic gas height zi/2 as a function of galactocentric radius after 
0.5 Gyr of evolution. Second row: fraction of gas beyond the characteristic height. Third row: vertical velocity dispersion. 
Solid lines: runs with H2. Dashed lines: runs without H2 



by the SPH kernel, so that the kicks are larger for particles 
closer to the new stellar particles. We have performed sim- 
ulations with various feedback efficiencies, helping to check 
this effect. Figure [20] shows the fraction of gas that is ex- 
pelled further than 1 kpc from the disc (or radially further 
than 60 kpc from the centre of the galaxy): without feed- 
back, the gas is only gravitationally heated and the effect of 
denser gas and higher clumping factor due to H2 makes the 
fraction of gas leaving the disc smaller than for a purely 
atomic hydrogen gas. However, with feedback, there is a 
higher fraction of gas outside the disc when H2 is included. 
This is because the higher concentration of the gas makes 
the feedback more efficient. In our simulations, the effect is 
especially striking in the outer parts of discs when H2 for- 
mation is taken into account, as the vertical restoring force 
is lower there. 

We checked that the vertical density distribution is con- 
sistent with both a more concentrated disc and a higher 
fraction of gas leaving the disc. Figure [2TI shows the verti- 
cal mass profile of the gas at large radii for = 40%: the 
gas is more concentrated in the disc plane, but the distribu- 
tion has higher density "tails" when more efficient feedback 
in denser regions allows the gas to be expelled from the 
disc. Figure 21 shows the characteristic height zi/2 of the 
gas, the distance from the disc plane for which the density 
equals half of the central density, as a function of radius 
for the various simulations. Especially at large radii, the 
characteristic height is lower for simulations with H2 for all 
feedbacks as H2 concentrates the gas in the middle plane. 
The fraction of gas beyond this height however increases 
with the feedback efficiency at large radii when H2 is taken 



into account, as gas is then efficiently expelled by feedback. 
Without feedback, some difference remains, probably due 
to gravitational heating produced by a higher clumping. We 
see indeed a slight difference in vertical velocity dispersion. 

3.2.5. Gas density profile 

The galaxies we have considered until now have a rather low 
gas surface density. We also run simulations with smaller 
characteristic radii, and therefore higher surface densities. 
Figure [22] displays the surface density for a Miyamoto- Nagai 
gas radius ig of 3.6 kpc. The transition radius between H2 
dominated and HI dominated regions, the radius at which 
Hh^ = Shi, has then a value similar to the average observed 
by Bigiel fc Blitz^ (2012} (their average value for nearby spi- 
ral galaxies observations is 14 Mq/pc^). Previous galaxies 
belong to the lower surface density group observed by BigielJ 
& Blitz (2012). The effect of H2 cooling is reduced in these 
galaxies with higher surface density, since a larger fraction 
of the gas belong to the central regions, with much more 
efficient star formation than the outer parts. H2 cooling is 
more important when there is more gas in the metal poor 
outer regions. 



4. Conclusions 

To explore the influence of molecular hydrogen in the 
physics of spiral discs, of star formation, and gas reservoirs 
in galaxy evolution, we have implemented in Gadget-2 de- 
tailed cooling by metals, for temperatures as low as 100 K, 
and cooling by H2 due to collisions with H, He and other H2 
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Fig. 22. Radial distribution of the surface density of H2, 
atomic gas HI and total hydrogen gas after 100 Myr for 
Qffb = 40% and a x x 500 UV scaling factor. 



molecules. The deter mination of the H2 density is inspired 
by the KMT recipe dKrumholz fc Gnedin|[2QTTl ) , usins; the 
stellar UV flux from young stars, and we study the influ- 
ence of cold and dense molecular phase together with stel- 
lar feedback on star formation. We have also implemented 
a stochastic star formation and a kinetic supernovae feed- 
back whose efficiency was varied in simulations including 
cooling by atomic/ionised gas and/or molecular hydrogen. 
The influence of H2 in the formation of dense gas and star 
formation is very important in the outer extended disk. 
Molecular hydrogen influences the vertical structure of the 
discs, especially when there is some stellar feedback: first 
H2 makes the gas more concentrated in the middle layer of 
the disc plane, but second the gas is also more susceptible 
of being ejected far from the disc, due to the higher effi- 
ciency of feedback in high density regions. Correlating SFR 
and gas surface density, it is found that molecular gas is a 
much better tracer of star formation than atomic g as, as 
is also observed in nearby galaxies. We find that including 
molecular hydrogen allows some slow star formation to oc- 
cur in the low metallicity outer parts of galaxies. If gas is 
accreted by the discs, it may help store some cold gas with 
a low star formation efficiency. 
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